3.物种分布数据清理

3.1 数据属性参数清理

3.1.1 使用is.na/ ! /dupicated

```{r eval=FALSE}

使用is.na和！同步排除

subset R&lonlat

rattler=subset(rattler, !is.na(lon) & !is.na(lat))

也可以使用is.na 来筛选错误记录；hh =subset(ratter,is.na(lon)) ; final_hh = ratter[!hh,]

查看筛选掉的数据量：

cat(nrow(occ_raw)-nrow(occ_clean), "records are removed")

使用duplicated去除重复

rattlerdups=duplicated(rattler[, c("lon", "lat")]) rattler <-rattler[!rattlerdups, ]

使用complete.case去除NAs值，通过此可以实现初步筛选，

occs <- occs[complete.cases(occs$longitude, occs$latitude), ]

另外一种清理分布数据方法或者说清理环境数据的方法：

Exlclude variables

cols.to.drop <- c('id', 'DWD_ID', 'STATION.NAME', 'FEDERAL.STATE', 'PERIOD') # columns to drop rows.to.drop <- complete.cases(dwd) # rows to drop dwd.data <- dwd[rows.to.drop, !(colnames(dwd) %in% cols.to.drop)]

#### 3.1.2 使用subset来清理数据集

```{r eval=FALSE} 
laurus <- gbif("Laurus", "nobilis")
##使用subset来筛选数据列
locs <- subset(laurus, select = c("country", "lat", "lon"))
head(locs)  # a simple data frame with coordinates
##subset条件筛选：
locs <- subset(locs, locs$lat < 90)
coordinates(locs) <- c("lon", "lat")  # set spatial coordinates
plot(locs)

3.1.3 基于达尔文核心来筛选

```{r eval=FALSE}

基于basisOfRecord包含多种观测数据类型，一般仅选择人类现今观测结果：

FOSSIL_SPECIMEN HUMAN_OBSERVATION(保留) LIVING_SPECIMEN
MACHINE_OBSERVATION OBSERVATION PRESERVED_SPECIMEN(保留) UNKNOWN

clean data

table(occ_unique$basisOfRecord) occ_unique_specimen <- subset(occ_unique, basisOfRecord=="PRESERVED_SPECIMEN")

这种打印方式比较有趣；

cat(nrow(occ_unique_specimen), "out of ", nrow(occ_unique), "records are specimen")

```{r eval=FALSE} 
##基于年代筛选数据
occ_final <- subset(occ_unique_specimen, year>=1950 & year <=2000)

3.1.4 核查物种分类单元

```{r eval=FALSE}

查看所下载的物种种类名称(条目)

sort(unique(occs$scientificName))

基于Taxonstand来核查所下载的物种数据是否属于正确分类单元，并且可以构成多物种数据框同步检查；

library(Taxonstand) species.names <- unique(occs$scientificName) tax.check <- TPL(species.names)

#### 3.1.5 核查物种的坐标系

```{r eval=FALSE} 
## 基于CoordinateCleaner包中提供的命令clean_coordinates
library(CoordinateCleaner)
## 第一步去除所有包含经纬度空值的数据行
occs.coord <- occs[!is.na(occs$decimalLatitude) & !is.na(occs$decimalLongitude),]
# 第二步使用clean_coordinates核查数据：仅输出正确数据；
geo.clean <- clean_coordinates(x = occs.coord,lon = "decimalLongitude",lat = "decimalLatitude",
species = "species",value = "clean")  ##value可调参
##第三步：可视化：
par(mfrow = c(1, 2))
plot(decimalLatitude ~ decimalLongitude, data = occs)
map(countriesLow,add = T)
plot(decimalLatitude ~ decimalLongitude, data = geo.clean)
map(countriesLow,add = T)

3.1.6 根据小数位数来去除低质量数据

## excle实现
=LEN(A1)-SEARCH(".",A1,1)
其中，len计算单元格的总长度，search用来查找小数点的位置，两者相减就是小数点的位数。

3.2 地理分布可视化筛选

3.2.1 分布点投影可视化

## 说明：
这一部分仅放置基本的wgs84格式的直接投影；如果需要更复杂的投影方式变换：参考栅格数据的投影；

```{r eval=FALSE}

R界面绘图可视化：

occ_final = read.csv("") coordinates(occ_final) <- ~ lon + lat ##这里的lon和lat是清理数据集中指定的经纬对应列 crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
proj4string(locs) <- crs.geo
summary(locs) data(wrld_simpl) plot(wrld_simpl) ## wrld_simpl是来自R内置的sp包； plot(locs, pch = 20, col = "steelblue")

```{r eval=FALSE} 
## R 界面web可视化：
## R::leaflets是一个用于人机交互的可视化地图包；
library(leaflet)
leaflet(s.sf.gcs) %>%  ## 其中s.sf.gcs是一个包含定义的栅格
  addPolygons() %>% 
  addTiles()

3.2.2 构筑范围可视化筛选(参见5.1.2)

```{r eval=FALSE}

说明：

构筑extent，除了如下代码所示利用范围筛选之外，最重要的一种方式即先确定边界范围（这通常可以对已知分布点求四象限求极值再加buffer确定）

extent的值是一个栅格分布的最大极值；具有顺序；

library(raster)

raster("bio2.tif") class : RasterLayer dimensions : 256, 360, 92160 (nrow, ncol, ncell) resolution : 0.04166667, 0.04166667 (x, y) extent : 107.841, 122.841, 20.16103, 30.8277 (xmin, xmax, ymin, ymax) crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 source : D:/zhang/bio2.tif names : bio2 ```

```{r eval=FALSE}

绘制范围上显示

data(wrld_simpl) data(worldHires)

构筑extend，注意这里是分布点极值—+1

plot(wrld_simpl, xlim=c(min(rattler$lon)-1,max(rattler$lon)+1), ylim=c(min(rattler$lat)-1,max(rattler$lat)+1), axes=TRUE, col="light yellow")

数据点图叠加显示

points(rattler$lon, rattler$lat, col="orange", pch=20, cex=0.75)

使用filter():直接筛选分割数据：

occs <- occs %>% filter(longitude < -20, longitude > -100, latitude < 12)

#### 3.2.3 使用google earth可视化点

```{r eval=FALSE}
##使用数据形式如下：
tips <- rep("xiz",47)
hh <- Andrew[,2:3]
LatLong <- paste0(hh[,1],":",hh[,2])
hh <- cbind(LatLong,tips)
hh <- as.data.frame(hh)
M1 <- gvisMap(hh, "LatLong", "tips", 
              options=list(showTip=TRUE, showLine=F, enableScrollWheel=TRUE, 
                           mapType='satellite', useMapTypeControl=TRUE, width=800,height=400))
plot(M1)

3.2.4 run(wallace)人机交互筛选

## 详见wallace## 详见(## 本书8.5章节wallace)

3.3 其他筛选与转换

3.3.3 计算点之间的距离

```{r eval=FALSE}

可以用于计算栅格的质心距离或者栅格质心转移

xy <- rbind(paris =c(89.5002419,29.27083354), c(90.14995087,28.95202496)) library(geosphere) ##注意投影坐标系的问题； distm(xy)

#### 3.3.4 shp转kml

```{r eval=FALSE} 
##其中以下的locs.gb需要包含完整的经纬度和坐标系
library(rgdal)
writeOGR(locs.gb, dsn = "locsgb.kml", layer = "locs.gb", driver = "KML")
##kml转shp
library(maptools)
newmap <- readOGR("locsgb.kml", layer = "locs.gb")
writePointsShape(locs.gb, "locsgb")
###栅格转kml
KML(tmin1.c, file = "tmin1.kml")

3.2.5 csv转shp

```{r eval=FALSE}

需要先把occ文件加上参考坐标投影，参见3.2.1

dir.create("temp") library(rgdal) shapefile(occ_final,"temp/occ_final.shp",overwrite=TRUE) loaded_shapefile <- shapefile("temp/occ_final.shp")

#### 3.3.6 百度坐标系转WGS84

```{r eval=FALSE} 
###百度转火星再转gs84；连续过程得到结果：
## 使用这个转化过程值需要定义工作空间和定义第一列的行名
## 然后在读取的csv文件中使用longitude和latitude
## 注意转换的最初数据格式已经转换为十进制；
setwd("C:/admin/desktop/")
data1 <- read.csv("three3.csv")
## head(data1) ##这里数据的列命名  data$longitude和data$latitude
##install.packages("data.table")
library(data.table)
x_pi <- 3.14159265358979324 * 3000.0 / 180.0
pi <- 3.1415926535897932384626#π
a <- 6378245.0#长半轴
ee <- 0.00669342162296594323#扁率

#百度坐标系(BD-09)转火星坐标系(GCJ-02)
data2<- data1[complete.cases(data1$longitude),]
num <- length(data2$longitude)  
library(data.table)  ##需加载data.table包
jwd_data <- data.table()
for(i in 1:num){
  data3 <- data1[i,]
  bd_lng <- data3$longitude
  bd_lat <- data3$latitude

  x <- bd_lng - 0.0065
  y <- bd_lat - 0.006
  z <- sqrt(x * x + y * y) - 0.00002 * sin(y * x_pi)
  theta <- atan2(y, x) - 0.000003 * cos(x * x_pi)
  lng <- z * cos(theta)
  lat <- z * sin(theta)

  #GCJ02(火星坐标系)转GPS84
  #纬度转换函数
  transformlat <- function(lng, lat) {
    ret <- -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat +
      0.1 * lng * lat + 0.2 * sqrt(abs(lng))
    + (20.0 * sin(6.0 * lng * pi) + 20.0 *
         sin(2.0 * lng * pi)) * 2.0 / 3.0
    + (20.0 * sin(lat * pi) + 40.0 *
         sin(lat / 3.0 * pi)) * 2.0 / 3.0
    + (160.0 * sin(lat / 12.0 * pi) + 320 *
         sin(lat * pi / 30.0)) * 2.0 / 3.0
    return(ret)
  }
  #纬经度转换函数
  transformlng <- function(lng, lat) {
    ret <- 300.0 + lng + 2.0 * lat + 0.1 * lng * lng +
      0.1 * lng * lat + 0.1 * sqrt(abs(lng))
    + (20.0 * sin(6.0 * lng * pi) + 20.0 *
         sin(2.0 * lng * pi)) * 2.0 / 3.0
    + (20.0 * sin(lng * pi) + 40.0 *
         sin(lng / 3.0 * pi)) * 2.0 / 3.0
    + (150.0 * sin(lng / 12.0 * pi) + 300.0 *
         sin(lng / 30.0 * pi)) * 2.0 / 3.0
    return(ret)
  }

  #火星坐标转换为WGS84
  dlat <- transformlat(lng - 105.0, lat - 35.0)
  dlng <- transformlng(lng - 105.0, lat - 35.0)
  radlat <- lat / 180.0 * pi
  magic <- sin(radlat)
  magic <- 1 - ee * magic * magic
  sqrtmagic <- sqrt(magic)
  dlat <- (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
  dlng <- (dlng * 180.0) / (a / sqrtmagic * cos(radlat) * pi)
  mglat <- lat + dlat
  mglng <- lng + dlng
  jwd_lng <- lng * 2 - mglng
  jwd_lat <- lat * 2 - mglat

  jwd_data1 <- data.table(jwd_lng,jwd_lat)
  jwd_data <- rbind(jwd_data,jwd_data1)
}
write.csv(jwd_data, "out_name.csv")  ##前面的jwd_data尽量不用改，后面导出数据命名可改；

3.3.7 CSV TO SHP

occ_final = read.csv("")
coordinates(occ_final) <- ~ lon + lat ##这里的lon和lat是清理数据集中指定的经纬对应列
crs.geo <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")  
proj4string(locs) <- crs.geo

3.3.8 SpatialPolygons to SpatialPolygonsDataFrame

## 其中range_py由分布范围直接生成：
occs <- na.exclude(py[,2:3])
colnames(occs) <- c("longitude","latitude")
xmin <- min(occs$longitude)-2
xmax <- max(occs$longitude)+2
ymin <- min(occs$latitude)-1
ymax <- max(occs$latitude)+3
bb <- matrix(c(xmin, xmin, xmax, xmax, xmin, ymin, ymax, ymax, ymin, ymin), ncol=2)
range_py <- sp::SpatialPolygons(list(sp::Polygons(list(sp::Polygon(bb)), 1)))

## 需要将空间数据框转为datafrme形式：
bc <- range_py
############################################
df <- data.frame(ID=character(), stringsAsFactors=FALSE )
for (i in bc@polygons ) { df <- rbind(df, data.frame(ID=i@ID, stringsAsFactors=FALSE))  }
# and set rowname=ID
row.names(df) <- df$ID
##############################################
rangepy <- SpatialPolygonsDataFrame(bc, df)
plot(rangepy)
## 将数据写出：

3.3.9 构建研究范围：

library(maptools)
data("wrld_simpl")
ranna <- frange(xh_na) 
poly_ployDF <- function(x){
  bc <- x
  df <- data.frame(ID=character(), stringsAsFactors=FALSE )
  for (i in bc@polygons ) { df <- rbind(df, data.frame(ID=i@ID, stringsAsFactors=FALSE))  }
  # and set rowname=ID
  row.names(df) <- df$ID
  rangepy <- SpatialPolygonsDataFrame(bc, df)
  return(rangepy)
  plot(rangepy)
}
ranna2<- poly_ployDF(ranna)
polygone1 <- rgeos::gBuffer(ranna2 , byid=TRUE, width=0)
polygone2 <- rgeos::gBuffer(wrld_simpl, byid=TRUE, width=0)
clip2 <- rgeos::gIntersection(polygone1,polygone2, byid=TRUE, drop_lower_td = TRUE)

3.4 物种分布数据降低采样偏差

```{r eval=FALSE}

使用spThin基于地理空间范围来筛选空间栅格内的数据点；

单位是千米；

library(spThin)

latitude/longitude must be in the WGS1984 pseudoprojection!

input <- read.csv("freckled madtom.csv")

thin to no more than one sample per kilometer##单位是千米！

If you're using on your own data, change 'out.dir', 'out.base', and 'log.file' to be what

change 'lat.col', 'long.col', and 'spec.col' to match your data

thin(loc.data = input, lat.col = "Lat", long.col = "Long", spec.col = "Common.Name",
thin.par = 1, reps = 10, locs.thinned.list.return = TRUE, write.files = TRUE, max.files = 1, ##max.file=1即可，输出一个用于建模 out.dir = "noct/", out.base = "noct", ##输出文件名及文件路径； write.log.file = TRUE,log.file = "noct.txt" )##输出日志

环境空间中降低采样偏差：

后两个步骤是使用pbdb R包的EnvSample函数执行的（Varela等人，2014年）。

### 3.5 物种分布数据清理相关R包

```r
## 由saramortara/data_cleaning 提供的基本数据清理技巧：
https://github.com/saramortara/data_cleaning

## 除一下专门开发用于清理物种分布的R包之外，在workflow阶段的wallace包也同样提供优秀的数据清理能力；

3.5 .1R:bdverse(shiny)

## bdverse基于R的shiny应用，用于批量清理物种分布数据，以及提供下载数据平台；
## 部分数据下载平台可能需要账号注册，以生成引用；
## bdverse包含三个部分：
1、bdDwc：提供交互式Rshiny界面；按照达尔文核心标准清理和下载数据
2、bdchecks:是用于执行，过滤，创建和管理各种生物多样性数据检查的基础架构;
3、bdclean：面向缺乏经验的R用户的用户友好型数据清理Shiny应用程序。它提供了用于管理完整的生物多样性数据清洗工作流程的功能，从收集用户输入以调整清洗程序到生成各种数据的报告和版本

## git网址：
https://bd-r.github.io/The-bdverse/
## R:bdverse使用教程：
https://bd-r.github.io/bdverse-user-guide/bddwc.html#bddwc-darwin-core
## R包使用：
# 整体启动：(##容安装出现问题)
install.packages("remotes")
remotes::install_github("bd-R/bdverse")
bdverse::bd_launcher()  ## 此种方式无法加载；

## 运行bdverse
bdverse_app()

3.5.2 mopa(scripts)

## mopa包的主要工具用途：
数据格式和预处理
伪缺席数据生成
模型拟合和预测
评估模型性能的其他选择
SDM预测分析

## R包安装：
if (!require(devtools)) install.packages("devtools")
devtools::install_git("https://github.com/SantanderMetGroup/mopa.git")
## R包的工作文档：
https://github.com/SantanderMetGroup/mopa/wiki

3.5.3 rangeModelMetadata

## rangeModelMetadata 的主要用途:
## 主要提供一系列数据清理工作：提供流程清理和多物种清理流程：
## R包安装：
devtools::install_github('cmerow / rangeModelMetaData')
## R包教程：
https://github.com/cmerow/rangeModelMetadata/tree/master/vignettes